rm(list=ls(all=TRUE))
setwd("~/Dropbox/Politics of Supreme Court Nominations/Chapter_Understanding_Yeas_and_Nays/PRSM_final_submission_and_replication_files")

options(digits=3)
library(tidyverse)
library(ggrepel)
library(readstata13)
library(sensemakr)
library(haven)


########################################################################################
#function for loess of binary data
binomial_smooth <- function(...) {
  geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
#GRAPH PARAMETERS (to use as shortcuts)
title.size <- 20
nom.label.size <- 10
x.axis.tick.size <- 14
y.axis.tick.size <- 14
x.axis.label.size <- 18
y.axis.label.size <- 18
strip.text.size <- 14
line.width <- 1.5
axis.text <- 1.2
label.text <- 1.2
y.offset <- .5
text.size <- 1.4
text.size2 <- 1.2
arrow.width <- 1.5
width <-2
text.size <- .7
y.axis.max <- 70
strip.color <- "black"
strip.text.size <- 20
########################################################################################

#READ IN DATASETS
#dataset of every nominee (Supreme Court Compendium) [ do not worry about factor label warnings; not relevant]
nominees <- read.dta13("sc_compendium_data_updated_all_nominees.dta", convert.underscore=TRUE, generate.factors=FALSE)
#newspaper coverage data
newspapers <- read.dta13("newspaper_data.dta")
#interest group participation data
interest_groups <- read.dta13("interest_group_data.dta")
#THOMAS DATA
#wide data
thomas.wide <- read.dta13("thomas_wide_data.dta", convert.underscore=TRUE)
#long data
thomas.long <- read.dta13( "thomas_long_data.dta", convert.underscore=TRUE)
#SOTOMAYOR AND KAGAN CCES DATA
cces.soto.kagan.wide <- read.dta13 ("CCES_soto_kagan_wide.dta", convert.underscore = TRUE) 
cces.soto.kagan.long <- read.dta13 ("CCES_soto_kagan_long.dta", convert.underscore = TRUE)

########################################################################################
#FIGURE 1: Percent yea votes over time
pdf("Figure_01.pdf", height=6, width=11)

ggplot( subset(nominees, nom.name!="J. Roberts (AJ)"), #drop J. Roberts 1st nomination as AJ
       aes(datenom,percent.yea.votes)) + theme_bw() + 
  geom_smooth() +   geom_hline(yintercept=.5, lty=2, col="yellow", lwd=1) +   ylab("Proportion of yea votes over time") +
  geom_point(data=subset(nominees, confirmed==0), aes(datenom,percent.yea.votes), shape=23, size=2, color="red") +   #do separate colors for confirmed and not
  geom_point(data=subset(nominees, confirmed==1), aes(datenom,percent.yea.votes), shape=19, size=2, color="dark green") + 
  geom_text_repel(data=subset(nominees, percent.yea.votes < 1), aes(datenom,percent.yea.votes, label = nom.name), size=5) +# add labels just for those with any nay votes
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),   
        axis.title.y = element_text(size=y.axis.label.size))  +
  scale_x_date(breaks=seq( as.Date('1/01/1800',format='%m/%d/%Y'), as.Date('1/01/2020',format='%m/%d/%Y'), "20 years"), date_labels =  "%Y") +
  scale_y_continuous(breaks=seq(0,1,.1)) +  coord_cartesian( ylim=c(.25,1)) 

dev.off()

########################################################################################
#Figure 2: Newspaper coverage of Supreme Court nominees, 1930-2020

pdf("Figure_02.pdf", height=6, width=11)

newspapers$date_announced <- as.Date(newspapers$date_announced, origin = "1960-01-01")#convert Stata dates to R

ggplot(newspapers , aes(y=n, date_announced)) +  theme_bw() +
  geom_point () +   geom_smooth(method='loess') +   labs(y="Number of newspaper stories", x="", title="") +
  geom_text_repel(data=newspapers,  aes(date_announced,n, label = newspapers$nom_name), size=3.5)+
  scale_x_date( breaks =  c( seq(as.Date("1930/1/1"), as.Date("2020/1/1"), "10 years")), labels=seq(1930,2020,10)) + 
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),   
        axis.title.y = element_text(size=y.axis.label.size))

dev.off()

########################################################################################
#Figure 3: Interest group mobilization on Supreme Court nominees, 1930-2020. 
pdf("Figure_03.pdf", height=6, width=11)

interest_groups$date_announced <- as.Date(interest_groups$date_announced, origin = "1960-01-01")#convert Stata dates to R

ggplot(interest_groups, aes(y=y, date_announced)) +  theme_bw() + 
  geom_point () +  geom_smooth(method='loess') +
  labs(y="Number of mobilizing interest groups", x="", title="") +
  geom_text_repel(data=interest_groups,  aes(date_announced,y, label = nom_name), size=3.5) +
  scale_x_date( breaks =  c( seq(as.Date("1930/1/1"), as.Date("2020/1/1"), "10 years")), labels=seq(1930,2020,10)) + 
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=14),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=14),   
        axis.title.y = element_text(size=18))

dev.off()

###########################################################################
###########################################################################
#RECALL ANALYSIS--LEVELS OF RECALL
###########################################################################
###########################################################################
#TABLE 1 (top portion)
  #THOMAS
100*prop.table(table(thomas.wide$index.recall.correct.thomas))#including DKs
100*prop.table(table(thomas.wide$index.recall.correct.thomas2)) #not including DKs
  #SOTO/KAGAN
100*prop.table(table(cces.soto.kagan.wide$index.recall.correct.nominee))#including DKs
100*prop.table(table(cces.soto.kagan.wide$index.recall.correct.nominee2))#not including DKs
#Table 1  (bottom portion)--correct responses  by party-line voting
  #THOMAS
thomas.long <- mutate(thomas.long,
                      sen.voted.with.party = ifelse( ( sen.dem==0 &  sen.vote==1 ) | #Republicans who voted for Thomas
                                                       ( sen.dem==1 &  sen.vote==0 ),1,0)   #Dems who voted against Thomas
)
#now compare recall to whether senator voted for party -- see need to account for NAs in individual resopoins
#merge index2 onto long data
foo <- dplyr::select(thomas.wide,case.id, index.recall.correct.thomas2)
thomas.long <- inner_join(thomas.long, foo)
thomas.party <- subset(thomas.long, !is.na(index.recall.correct.thomas2))#drop people who guessed for either nominee
#how often are respondents correct when senator votes against party?
with(thomas.party, 100*mean(sen.recall.correct.thomas, na.rm=TRUE) )#overall
with(thomas.party, 100*mean(sen.recall.correct.thomas[sen.voted.with.party==1], na.rm=TRUE) )
with(thomas.party, 100*mean(sen.recall.correct.thomas[sen.voted.with.party==0], na.rm=TRUE) )

#SOTO/KAGAN
#for RR, break down recall by whether senator voted against party line (use long data)
cces.soto.kagan.party <- mutate(cces.soto.kagan.long,
                                sen.voted.with.party = ifelse( ( sen.dem==0 &  sen.confirm.nominee==0 ) | #Republicans who against Soto/Kagan
           
                                                                                                                      ( sen.dem==1 &  sen.confirm.nominee==1),1,0)   #Dems who voted for Soto/Kagan
)
foo <- dplyr::select(cces.soto.kagan.wide,case.id.new, index.recall.correct.nominee2)
cces.soto.kagan.party <- inner_join(cces.soto.kagan.party, foo, by="case.id.new")
cces.soto.kagan.party <- subset(cces.soto.kagan.party, !is.na(index.recall.correct.nominee2))#drop people who guessed for either nominee
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee, na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[sen.voted.with.party==1], na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[sen.voted.with.party==0], na.rm=TRUE) )

####################################################################################
#Table A-5: Additional correct recall analysis
####################################################################################
#)Voters are flipping weighted coins with weights given by the actual number of yes vs no votes on a justice. 
#thomas was confirmed 52-48
set.seed(212323)
thomas.wide <- thomas.wide %>% mutate(
  guess.vote1 = rbinom(nrow(thomas.wide), 1, .52))
set.seed(33)
thomas.wide <- thomas.wide %>% mutate(
  guess.vote2 = rbinom(nrow(thomas.wide), 1, .52),
  guess1.correct = ifelse(  guess.vote1 == sen1.vote,1,0),
  guess2.correct = ifelse(  guess.vote2 == sen2.vote,1,0),
  guess.index = guess1.correct + guess2.correct
)
#COUNTERFACTUAL DISTRIBUTION [actual reported above]
100*prop.table(table(thomas.wide$guess.index))

#sotomayor was confirmed 68-31 /kagan 63-37
soto.prop <- (68/99)
kagan.prop <- .67
foo <- subset( cces.soto.kagan.wide, !is.na(sen2.confirm.nominee))#drop obs where sen2.confirm.nominee has NA (FL and Mass due to senator turnover/absentia)
set.seed(214)
foo <- foo %>% mutate(
  guess.vote1 = case_when(
    nominee=="sotomayor" ~  rbinom(nrow(foo), 1, soto.prop),
    nominee=="kagan" ~  rbinom(nrow(foo), 1, kagan.prop)))
set.seed(319)
foo<- foo%>% mutate(
  guess.vote2 = case_when(
    nominee=="sotomayor" ~  rbinom(nrow(foo), 1, soto.prop),
    nominee=="kagan" ~  rbinom(nrow(foo), 1, kagan.prop)),
  guess1.correct = ifelse(  guess.vote1 == sen1.confirm.nominee,1,0),
  guess2.correct = ifelse(  guess.vote2 == sen2.confirm.nominee,1,0),
  guess.index = guess1.correct + guess2.correct
)
#COUNTERFACTUAL DISTRIBUTION [actual reported above]
100*prop.table(table(foo$guess.index))#distribution under guessing. [actual reported above]

##Table A-17--separate results for Sotomayor and Kagan
  #top-portion
100*prop.table(table(cces.soto.kagan.wide$index.recall.correct.nominee[cces.soto.kagan.wide$nominee=="sotomayor"]))#including DKs
100*prop.table(table(cces.soto.kagan.wide$index.recall.correct.nominee2[cces.soto.kagan.wide$nominee=="sotomayor"]))#not including DKs
100*prop.table(table(cces.soto.kagan.wide$index.recall.correct.nominee[cces.soto.kagan.wide$nominee=="kagan"]))#including DKs
100*prop.table(table(cces.soto.kagan.wide$index.recall.correct.nominee2[cces.soto.kagan.wide$nominee=="kagan"]))#not including DKs
  #bottom portion
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[nominee=="sotomayor"], na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[sen.voted.with.party==1 & nominee=="sotomayor"], na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[sen.voted.with.party==0 & nominee=="sotomayor"], na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[nominee=="kagan"], na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[sen.voted.with.party==1 & nominee=="kagan"], na.rm=TRUE) )
with(cces.soto.kagan.party, 100*mean(sen.recall.correct.nominee[sen.voted.with.party==0 & nominee=="kagan"], na.rm=TRUE) )


###########################################################################
###########################################################################
#VOTER RECALL AND POLITICAL ENGAGEMENT: figures 4 and 5
###########################################################################
###########################################################################

#Figure 4: Political engagement and voter recall of senators’ votes to confirm Clarence Thomas.
thomas.wide <- mutate(thomas.wide,  partisanship.factor = ifelse(dem.lean==1, "2_DEM.", ifelse( gop.lean==1, "3_REP.", "1_IND")),  partisanship.binary = factor(ifelse(dem.lean == 1 | gop.lean ==1,"Partisan","Independent")))
thomas.wide.cut <- dplyr::select(thomas.wide, index.recall.correct.thomas,  pol.knowledge,  ideo.extremity.7, pol.attn,  educ,  partisanship.binary, media.index.quart, R.voted.gen.elec)
foo.wide <- gather(thomas.wide.cut, variable, level, -index.recall.correct.thomas)
head(foo.wide)
foo.thomas.recall.means <- foo.wide %>% group_by(variable, level) %>% summarise( mean.recall = mean(index.recall.correct.thomas))
foo.thomas.recall.means = filter(foo.thomas.recall.means, !is.na(level))
foo.thomas.recall.means$variable <- factor(foo.thomas.recall.means$variable)
levels(foo.thomas.recall.means$variable) <- c("Education", "Ideological extremity", "Media consumption", "Partisan identification", "Political attention", "Political knowledge", "Voted in presidential election")

pdf("Figure_04.pdf", height=11, width=12)

ggplot(foo.thomas.recall.means, aes(level,mean.recall)) + theme_bw() + 
  geom_bar( position = "dodge", stat="identity", width = 0.2) +  facet_wrap(~variable, scales="free_x", ncol=2) + 
  xlab("Level of each predictor, from low engagement to high") +  ggtitle("Thomas") + 
  ylab("Mean index of Thomas recall") + #scale_y_continuous(expand = c(0, 0)) +
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(angle = 0, hjust = .5,size=x.axis.tick.size),axis.title.x = element_text(size=x.axis.label.size),
        axis.text.y = element_text(size=y.axis.tick.size),  axis.title.y = element_text(size=y.axis.label.size),
                strip.text.x = element_text(size = strip.text.size, colour = strip.color))

dev.off()

#Figure 5,  Political engagement and voter recall of senators’ votes to confirm Sotomayor/Kagan
#use vars related to voter engagement
cces.soto.kagan.wide.cut <- dplyr::select(cces.soto.kagan.wide, nominee, index.recall.correct.nominee,news.interest,ideo.extremity.7,educ, partisanship.binary,  knowledge.index.party.control, R.voted.gen.elec, registered.to.vote)
#drop missing responses to recall
cces.soto.kagan.wide.cut <- filter(cces.soto.kagan.wide.cut ,!is.na( index.recall.correct.nominee))
cces.soto.kagan.wide.cut$partisanship.binary <- ifelse( cces.soto.kagan.wide.cut$partisanship.binary =="", NA, cces.soto.kagan.wide.cut$partisanship.binary)
foo.gather <- gather(cces.soto.kagan.wide.cut, variable, level, -index.recall.correct.nominee, -nominee)
#create means by group, first pooling soto and kagan
soto.kagan.pooled.engage <- foo.gather %>% group_by(variable, level) %>% summarise(
  mean.recall = mean(index.recall.correct.nominee))
soto.kagan.pooled.engage = filter(soto.kagan.pooled.engage, !is.na(level))
soto.kagan.pooled.engage$variable <- factor(soto.kagan.pooled.engage$variable)
levels(soto.kagan.pooled.engage$variable) <- c("Education", "Ideological extremity", "Political knowledge", "News interest", "Partisan identification", "Voted in presidential election", "Registered to vote")

pdf("Figure_05.pdf", height=11, width=12)

ggplot( soto.kagan.pooled.engage, aes(level,mean.recall)) +   theme_bw() + 
  geom_bar( position = "dodge", stat="identity", width = 0.1) + facet_wrap(~variable, scales="free_x", ncol=2) + xlab("Level of each predictor, from low engagement to high") +
  ggtitle("Sotomayor and Kagan") + 
  ylab("Mean index of nominee recall (Sotomayor and Kagan pooled)") + #scale_y_continuous(expand = c(0, 0)) +
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(angle = 0, hjust =.51,size=x.axis.tick.size),axis.title.x = element_text(size=x.axis.label.size),
        axis.text.y = element_text(size=y.axis.tick.size),  axis.title.y = element_text(size=y.axis.label.size),
        strip.text.x = element_text(size = strip.text.size, colour = strip.color))

dev.off()

########################################################################
#Table 5: Approval of senators by PERCEIVED nominee and party agreement. 
########################################################################
#THOMAS
thomas.mean.approval.by.perceptions <- thomas.long %>% group_by ( sen.perceived.party.agreement, sen.perceived.thomas.agreement) %>%   summarise(
  mean.approval = mean(sen.approval.binary,na.rm=TRUE))
arrange(thomas.mean.approval.by.perceptions, -sen.perceived.party.agreement, -sen.perceived.thomas.agreement)

#calculate marginals and overall mean
thomas.long %>%  group_by (sen.perceived.thomas.agreement) %>%   summarise(
  mean.approval = mean(sen.approval.binary,na.rm=TRUE))
thomas.long %>%  group_by (sen.perceived.party.agreement) %>%   summarise(
  mean.approval = mean(sen.approval.binary,na.rm=TRUE))
mean(thomas.long$sen.approval.binary,na.rm=TRUE) #overall

#SOTO/KAGAN
sk.mean.approval.by.perceptions <- cces.soto.kagan.long %>% group_by ( sen.perc.party.agreement, sen.perc.nominee.agreement) %>%   summarise(
  mean.approval = mean(sen.approval.binary,na.rm=TRUE))
sk.mean.approval.by.perceptions <- filter(sk.mean.approval.by.perceptions , !is.na(sen.perc.nominee.agreement) & !is.na(sen.perc.party.agreement ))
arrange(sk.mean.approval.by.perceptions , -sen.perc.party.agreement, -sen.perc.nominee.agreement)

#calculate marginals and overall mean
cces.soto.kagan.long %>%  group_by (sen.perc.nominee.agreement) %>%   summarise(
  mean.approval = mean(sen.approval.binary,na.rm=TRUE))
cces.soto.kagan.long %>%  group_by (sen.perc.party.agreement) %>%   summarise(
  mean.approval = mean(sen.approval.binary,na.rm=TRUE))
mean(cces.soto.kagan.long$sen.approval.binary,na.rm=TRUE)#overall 

############################################################################################
#MAKE FIGURE 7
  #note the do file must first be run to create "saved_regression_results_for_figure7.dta for this code to work
############################################################################################
regs <- read.dta13("saved_regression_results_for_figure7.dta")
table(regs$label)

#APPROVAL
regs.approval <-  regs %>%  filter(str_detect(label, 'iv_with_controls') & str_detect(label, 'approval') & str_detect(var, 'sen_perc_') &  var!="sen_perc_party_agreement")
table(regs.approval$label)
regs.approval <- arrange(regs.approval, var)
regs.approval 

regs.approval$graph.label <- c("ACA", "CHIP", "Dodd-Frank", "Don't Ask, Don't Tell",  "Lilly Ledbetter", "Sotomayor/Kagan", "Stimulus")
regs.approval <-  dplyr::select(regs.approval, var, coef,stderr, graph.label)
regs.approval <-arrange(regs.approval, -coef)
y.approval <- c(nrow(regs.approval):1)
#VOTE CHOICE
regs.vote.choice <-  regs %>%  filter(str_detect(label, 'iv_with_controls') & str_detect(label, 'vote_choice') & str_detect(var, 'sen_perc_') &  var!="sen_perc_party_agreement") 
table(regs.vote.choice$label)
regs.vote.choice <- arrange(regs.vote.choice, var)
regs.vote.choice 

regs.vote.choice$graph.label <- c("ACA","Dodd-Frank", "Don't Ask, Don't Tell", "Kagan", "Stimulus")
regs.vote.choice <-  dplyr::select(regs.vote.choice, var, coef,stderr, graph.label)
regs.vote.choice <-arrange(regs.vote.choice, -coef)
y.vote.choice <- c(nrow(regs.vote.choice):1)

pdf("Figure_07.pdf", height=11, width=9)

layout(c(1,2), heights=c(1.4,1))
par(mar=c(5,13,2,2))
plot(regs.approval$coef,y.approval, xlim=c(0,.3), type="n", axes=FALSE, xlab="", ylab="", xaxs="i")
abline(h=y.approval, lty=2, col="gray30")
segments(regs.approval$coef-1.96*regs.approval$stderr,y.approval,  regs.approval$coef+1.96*regs.approval$stderr, y.approval, lwd=2)
points(regs.approval$coef,y.approval, pch=21, bg="white", cex=2)
axis(1, mgp=c(2,.5,0), cex.axis=1.5)
axis(2, at = y.approval, labels=regs.approval$graph.label, las=2,  cex.axis=1.5)
mtext("Effect of perceived agreement on senator approval", 1, line=2.3, cex=1.7, xpd=NA)
box()
#draw box around Soto/Kagan
x.end <- -.009
x.begin <- -.11
polygon(x=c(x.begin,x.begin, x.end, x.end), y=c(7.3, 6.7, 6.7, 7.3), xpd=NA, border="red", lwd=2)
mtext("Approval",3, cex=2, font=2, line=.5)

plot(regs.vote.choice$coef,y.vote.choice, xlim=c(0,.3), type="n", axes=FALSE, xlab="", ylab="", xaxs="i")
abline(h=y.vote.choice, lty=2, col="gray30")
segments(regs.vote.choice$coef-1.96*regs.vote.choice$stderr,y.vote.choice,  regs.vote.choice$coef+1.96*regs.vote.choice$stderr, y.vote.choice, lwd=2)
points(regs.vote.choice$coef,y.vote.choice, pch=21, bg="white", cex=2)
axis(1, mgp=c(2,.5,0), cex.axis=1.5)
axis(2, at = y.vote.choice, labels=regs.vote.choice$graph.label, las=2,  cex.axis=1.5)
mtext("Effect of perceived agreement on senator vote choice", 1, line=2.3, cex=1.7, xpd=NA)
box()
#draw box around Kagan
x.end <- -.009
x.begin <- -.05
polygon(x=c(x.begin,x.begin, x.end, x.end), y=c(3.3, 2.7, 2.7, 3.3), xpd=NA, border="red", lwd=2)
mtext("Vote Choice",3, cex=2, font=2, line=.5)

dev.off()

############################################################################################
#COVARIATE BALANACE ANALYSIS IN APPENDIX (Table A-1)
############################################################################################
cces.soto.kagan.long <- mutate(cces.soto.kagan.long
                               , treated =  case_when( sen.actual.nominee.agreement == 1 ~ 1, sen.actual.nominee.agreement == -1 ~ 0 )
                               , educ.level1 = ifelse(educ==1,1,0)
                               , educ.level2 = ifelse(educ==2,1,0)   
                               , educ.level3 = ifelse(educ==3,1,0)   
                               , educ.level4 = ifelse(educ==4,1,0)
                               , age.level1 = ifelse(age.cat==1,1,0)
                               , age.level2 = ifelse(age.cat==2,1,0)   
                               , age.level3 = ifelse(age.cat==3,1,0)   
                               , age.level4 = ifelse(age.cat==4,1,0)
                               , white = ifelse(race.wbh==1,1,0)
                               , black = ifelse(race.wbh==2,1,0)
                               , hispanic = ifelse(race.wbh==3,1,0)
)
head(cces.soto.kagan.long)

soto.kagan.balance <- dplyr::select(cces.soto.kagan.long, sen.actual.nominee.agreement, dem.lean, ind.lean, gop.lean,  white, black,hispanic, contains("level"),engage.factor.resc )
#create means by agreeement
soto.kagan.means <- soto.kagan.balance %>% group_by(sen.actual.nominee.agreement) %>% 
  summarise_all ( ~mean(.x, na.rm = TRUE))
round(soto.kagan.means,2)
#reshape to long, with one row for each variable
soto.kagan.means  <- pivot_longer(soto.kagan.means, - sen.actual.nominee.agreement, names_to="variable", values_to="mean")
soto.kagan.means <- arrange(soto.kagan.means, variable)
soto.kagan.means 
soto.kagan.means <- as.data.frame( spread(soto.kagan.means,   sen.actual.nominee.agreement,  mean))
#soto.kagan.means is what is reported in table means
soto.kagan.means$foo <- c(1,2,3,4, 15, 5, 13, 7,8,9,10, 11, 6, 14, 12)
soto.kagan.means <- soto.kagan.means[ordered(soto.kagan.means$foo),]
soto.kagan.means$foo<-NULL
options(digits=2)
soto.kagan.means 

############################################################################################
#SENSITIVITY ANALYSIS
############################################################################################
#create new version of long data with simplified names for party and ideological agreement
cces.data <- read_dta("CCES_soto_kagan_long.dta") %>% 
  rename(Copartisan =  sen_actual_party_agreement,  Coideology = sen_actual_ideo_agreement )

#vote choice
cces.data.vc <- subset(cces.data, nominee =="kagan" & R_had_incumbent10_elec==1) 
#need to drop incumbents who didn't run or lost in primary
cces.data.vc <- subset(cces.data.vc, ! sen_cces_name  %in%      
                         c( "Christopher Dodd" ,"George LeMieux" ,"Roland Burris" ,"Evan Bayh" ,"Samuel Brownback" ,
                            "Jim Bunning" ,"Christopher 'Kit' Bond" ,"Judd Gregg" ,"Byron Dorgan" ,"George Voinovich"  
                            ,"Arlen Specter" ,"Robert Bennett"))
#FIRST STAGE
#APPROVAL
fit_cces_fs_approval <- lm( sen_perc_nominee_agreement ~ sen_actual_nominee_agreement + Copartisan 
                            + Coideology  +  as.factor(educ) + female + as.factor(race_wbh) + as.factor(age_cat) + dem_lean + gop_lean  + engage_factor, 
                            data=cces.data)
arm::display(fit_cces_fs_approval, detail=TRUE)

sens_cces_fs_approval_party <- sensemakr(fit_cces_fs_approval, 
                                         treatment = "sen_actual_nominee_agreement",
                                         benchmark_covariates = "Copartisan",
                                         kd = c(1,2))
summary(sens_cces_fs_approval_party, digits=2)

plot( sens_cces_fs_approval_party , col.contour = "darkgray", frame =TRUE, nlevels = 10,
      main = "First stage, Approval/Copartisan",    label.bump.x = 0.03,    cex.label.text = 0.75)

sens_cces_fs_approval_ideology <- sensemakr(fit_cces_fs_approval, 
                                            treatment = "sen_actual_nominee_agreement",
                                            benchmark_covariates = "Coideology",
                                            kd = c(5,10))
summary(sens_cces_fs_approval_ideology, digits=2)

plot( sens_cces_fs_approval_ideology  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "First stage, Approval/Ideology",    label.bump.x = 0.03,    cex.label.text = 0.75)

#VOTE CHOICE
fit_cces_fs_vc <- lm( sen_perc_nominee_agreement ~ sen_actual_nominee_agreement + Copartisan 
                      + Coideology  +  as.factor(educ) + female + as.factor(race_wbh) + as.factor(age_cat) + dem_lean + gop_lean  + engage_factor, 
                      data=cces.data.vc)
arm::display(fit_cces_fs_vc, detail=TRUE)

sens_cces_fs_vc_party <- sensemakr(fit_cces_fs_vc, 
                                   treatment = "sen_actual_nominee_agreement",
                                   benchmark_covariates = "Copartisan",
                                   kd = c(1,2))
summary(sens_cces_fs_vc_party, digits=2)

plot( sens_cces_fs_vc_party , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "First stage, Vote Choice/Copartisan",    label.bump.x = 0.03,    cex.label.text = 0.75)

sens_cces_fs_vc_ideology <- sensemakr(fit_cces_fs_vc, 
                                      treatment = "sen_actual_nominee_agreement",
                                      benchmark_covariates = "Coideology",
                                      kd = c(5,10))
summary(sens_cces_fs_vc_ideology, digits=2)

plot( sens_cces_fs_vc_ideology  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "First stage, Vote Choice/Ideology",    label.bump.x = 0.03,    cex.label.text = 0.75)

#make first-stage plot (Figure A_1)
pdf("Figure_A1.pdf", height=8, width=8)

par(mfrow=c(2,2), mar=c(2,4,2,2))

plot( sens_cces_fs_approval_party , col.contour = "darkgray", frame =TRUE, nlevels = 10,
      main = "First stage, Approval/Copartisan",    label.bump.x = 0.03,    cex.label.text = 0.75)
plot( sens_cces_fs_approval_ideology  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "First stage, Approval/Ideology",    label.bump.x = 0.03,    cex.label.text = 0.75)
plot( sens_cces_fs_vc_party , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "First stage, Vote Choice/Copartisan",    label.bump.x = 0.03,    cex.label.text = 0.75)
plot( sens_cces_fs_vc_ideology  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "First stage, Vote Choice/Ideology",    label.bump.x = 0.03,    cex.label.text = 0.75)

dev.off()

######################################################################
#Reduced form
######################################################################
#Approval
#Partisanship
fit_cces_rf_approval <- lm(sen_approval_binary~ sen_actual_nominee_agreement +  Copartisan + Coideology +
                             + Coideology +  as.factor(educ) + female + as.factor(race_wbh) + as.factor(age_cat) + dem_lean + gop_lean  + engage_factor, 
                           data=cces.data)
arm::display(fit_cces_rf_approval, detail=TRUE)

sens_cces_rf_approval_party <- sensemakr(fit_cces_rf_approval,treatment = "sen_actual_nominee_agreement", benchmark_covariates = "Copartisan",
                                         kd = c(.5,1,2))
summary( sens_cces_rf_approval_party , digits=2)

plot( sens_cces_rf_approval_party  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression (Approval/Copartisan)",    label.bump.x = 0.03,    cex.label.text = 0.75)

#Ideology
sens_cces_rf_approval_ideology <- sensemakr(fit_cces_rf_approval,
                                            treatment = "sen_actual_nominee_agreement", benchmark_covariates = "Coideology",
                                            kd = c(1,5))
summary(sens_cces_rf_approval_ideology , digits=2)

plot( sens_cces_rf_approval_ideology  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression (Approval/Coideology)",    label.bump.x = 0.03,    cex.label.text = 0.75)
#Vote Choice
#partisanship
fit_cces_rf_vote <- lm(R_vote_inc ~ sen_actual_nominee_agreement +  Copartisan + Coideology +
                         + Coideology +  as.factor(educ) + female + as.factor(race_wbh) + as.factor(age_cat) + dem_lean + gop_lean  + engage_factor, 
                       data=cces.data.vc)
arm::display(fit_cces_rf_vote, detail=TRUE)

sens_cces_rf_vote_party <- sensemakr(fit_cces_rf_vote, 
                                     treatment = "sen_actual_nominee_agreement",
                                     benchmark_covariates = "Copartisan",
                                     kd = c(.5))
summary(sens_cces_rf_vote_party, digits=2)

plot( sens_cces_rf_vote_party , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression (Vote choice)",    label.bump.x = 0.03,    cex.label.text = 0.75)

sens_cces_rf_vote_ideology <- sensemakr(fit_cces_rf_vote, 
                                        treatment = "sen_actual_nominee_agreement",
                                        benchmark_covariates = "Coideology",
                                        kd = c(1,5))
summary(sens_cces_rf_vote_ideology, digits=2)

plot( sens_cces_rf_vote_ideology , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression (Vote choice)",    label.bump.x = 0.03,    cex.label.text = 0.75)

###################################################
#Reduced form (Figure A-2)
###################################################
pdf("Figure_A2.pdf", height=8, width=8)

par(mfrow=c(2,2), mar=c(2,4,2,2))


plot( sens_cces_rf_approval_party  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression, (Approval/Copartisan)",    label.bump.x = 0.03,    cex.label.text = 0.75)

plot( sens_cces_rf_approval_ideology  , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression, (Approval/Coideology)",    label.bump.x = 0.03,    cex.label.text = 0.75)

plot( sens_cces_rf_vote_party , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression, (Vote choice/Copartisan)",    label.bump.x = 0.03,    cex.label.text = 0.75)

plot( sens_cces_rf_vote_ideology , col.contour = "darkgray", frame = TRUE, nlevels = 10,
      main = "Outcome regression, (Vote choice/Coideology)",    label.bump.x = 0.03,    cex.label.text = 0.75)

dev.off()


sink()















